function z = Du(p)
    x = p(:,1); y = p(:,2);
    z(:,1) = (2.*x + 1)./((x + 1/2).^2 + (y + 1/2).^2 + 1/100).^2 ...
                - (2.*x - 1)./((x - 1/2).^2 + (y - 1/2).^2 + 1/100).^2;
    z(:,2) = (2.*y + 1)./((x + 1/2).^2 + (y + 1/2).^2 + 1/100).^2 ...
                - (2.*y - 1)./((x - 1/2).^2 + (y - 1/2).^2 + 1/100).^2;
end
